home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
sound
/
fftscop4.zip
/
HILFE.TXT
< prev
next >
Wrap
Text File
|
1994-05-30
|
10KB
|
242 lines
SMALLFFT.EXE eine Abgemagerte Version um die Programmierung
SMALLFFT.C besser verstehen zu können!
FFT.ZIP die FFT-Routinen, alt (1988) aber gut!
jedoch mangelhaft dokumentiert...
Kurzbeschreibung FFT.ZIP
---------------------------------------------------------------------------
die aufgenommenen Werte werden im real[] Array abgelegt,
das imag[] Array wird mit Nullen aufgefüllt...
dann werden die Globalen Variablen samples und power gesetzt:
samples muß ein Potenz von 2 sein: 2,4,8,16,32,64,128,256,512,1024,2048,...
int samples, power;
double real[2048], imag[2048], max;
samples = Anzahl aufgenommener Werte im real[]/imag[] Array;
power = log10((double)samples) / log10((double)2.0);
Jetzt wird die Funktion fft() ausgeführt!
dieses ist sehr schnell, extrem viel schneller als die normale FT.
das Real/Imaginär-Spektrum befindet sich danach im real[]/imag[] Array.
Allerdings recht merkwürdig durcheinandergewürfelt (das ist bei FFT normal)!
die permute(n) liefert uns die Position des Spektrums n in real[]/imag[]
Beispiel:
t=permute(n)
r=real[t]
i=imag[t]
Betrag=sqrt(r*r + i*i)
return Betrag;
dieses Beispiel ist die Funktion magnitude(n), welche den Betrag liefert!
Wenn wir uns eh nur für das Betragspektrum interessieren benutzen wir
also magnitude() und haben mit permute() nix mehr am Hut!
Iss auch ne Rückwärts (inverse) Fast-Fourier-Transfo möglich?
Aber klaro, laut dem schlauen FFT-Transformation Buch das ich mal kurz
in den Fingern hatte ist dazu lediglich der Imaginärteil zu invertieren,
(also mit -1 zu multiplizieren), fft() aufrufen, und dann wieder den
Imaginärteil invertieren.
aber dabei ist die Richtige Reihenfolge im Array wichtig (siehe permute)!
fft();
for(i=0;i<samples;i++) imag[i]*=(-1);
fft();
for(i=0;i<samples;i++) imag[i]*=(-1);
nach diesem kurzen Programm stehen nicht etwa wieder die ursprünglichen
Werte im Array, sondern Datenmüll, weil die Verwürfelung ignoriert wurde!
Richtiger wäre:
fft();
for(i=0;i<samples;i++)
{
t=permute(i)
real2[i]=real[t];
imag2[i]=imag[t];
}
for(i=0;i<samples;i++) real[i]=real2[i], imag[i]=imag2[i];
for(i=0;i<samples) imag[i]*=(-1);
fft();
for(i=0;i<samples) imag[i]*=(-1);
for(i=0;i<samples;i++)
{
t=permute(i)
real2[i]=real[t];
imag2[i]=imag[t];
}
Nach der Ausführung dieses Programs stehen die urspünglichen Werte mit
kleinen/grossen Fehlern wieder im real[]/imag[] Array und im
Zwischenspeicher-Array real2[]/imag2[] finden wir das Real/Imaginärspektrum!
welches Format hat denn nu dasch Real/Imaginär/Betragsspektrum?
0 = Gleichanteil ....... samples-1 = samplerate
leider können wir das Array nur bis zur Hälfte (samples/2 - 1) nutzen weil
ein gewisser Nyquist das sogenannte Nyquist-Kriterium aufgestellt hat,
nach dem bei Abtastung maximal Frequenzen der halben Abtastrate abgetastet
werden können.
Erklärung mittels Systemtheorie: (bahh, würg, igittigitt)
Ich versuche es unmathematisch visual anschaulich zu erklären,
also Zeichblatt und Stift rauskramen und versuchen das Folgende zu zeichnen:
wir betrachten die Abtastung im Freqenzbereich:
die kontinuierliche Abtastung (Kammfunktion im Frequenz- und Zeitbereich)
mit Abstand samplerate (Hz) im Freqenzbereich und
mit Abstand 1/samplerate (s) im Zeitbereich
wird im
Freqenzbereich gefaltet mit dem Eingangssignalspektrum und im
Zeitbereich multipliziert mit der Eingangsignal.
Faltung ist normalerweise besch... zu zeichnen aber wegen den Diracstößen
und der Kammfunktion wirds viel einfacher, wir brauchen bloss das
Eingangssignalspektrum um jeden Diracstoß zu zeichnen, als wenn jeder
Diracstoß f=0Hz hätte, wenn wir das gezeichnet haben,
verstehen wir nicht nur das Nyquist-Kriterium,
wir wissen auch wie unser Spektrum zwischen samples/2 ... samples aussieht!
Das sind nähmlich nicht weiter als die negativen Frequenzen des um
Samplerate gefalteten Eingangssignalspektrums.
Also entspricht der linke Rand (0) und der rechte Rand (samples) des Arrays
(permute beachten) dem Gleichanteil (0 Hz) des Eingangsspektrums, da der
maximale adressierbare Wert im Array jedoch samples-1 ist finden wir den
Gleichanteil nur bei real[permute(0)],imag[permute(0)]...
Realteil:
das Signal von 0...samples/2 - 1 wird um samples/2 gespiegelt.
Imaginärteil:
das Signal von 0...samples/2 - 1 wird um samples/2 gespiegelt, und zusätzlich
invertiert (also das Vorzeichen gewechselt).
Betragspektrum: um das Betragsspektrum in den Zeitbereich zu transformieren
zerlegen wir es in eine gerade und eine ungerade Funktion,
der gerade Anteil kommt nach real[]
der ungerade Anteil kommt nach imag[]
dann wird Realteil und Imagiärteil wie angegeben gespiegelt!
für die inverse FFT ist natürlich noch die weiter oben erwähnte doppelte
Invertierung des Imaginärteils notwendig!!!
-----------------------------------------------------------------------------
Informationstragend ist dabei eigendlich nur der Bereich von 0...samples/2-1
entsprechend f=0 bis f=samplerate/2 also von 0 Hz bis zur halben Abtastfrequenz.
Wie meine Experimente mittels Soundkarte und Funktionsgenrator zeigten,
ist beim Realspektrum und beim Imaginärspektrum das Vorzeichen dauernd am
kippen, schade...
das Betragsspektrum wird mittels Arctan(imag[permute(i)]/real[permute(i)])
berechnet, leider eine ungerade und damit vorzeichensensitive Funktion...
das Betragsspektrum ist nicht nur sehr aussagekräftig, es ist auch wegen
der Quadratbildung nicht vorzeichensensitiv also stabil!!!
es gibt da allerdings noch ein paar Probleme auf die ich gestoßen bin:
bedingt durch die Abtastung kommt es wenn sich ein Vielfaches der
Abgetasteten Frequenz sich der Abtastrate nähert zu einer Schwebung
(Amplitudenmodulation), daraus resultieren mit Sicherheit Fehler!
die FFT hat auch noch ein ähnlich gelagertes Problem:
wie wir aus unserer Zeichnung ersehen können, ist die diskrete
Fourier Transformation bedingt durch die Abtastung
im Freqenzbereich periodisch, daraus ergibt sich der Umkehrschluß das
auch irgendeine Art von Periodischen Verhalten im Zeitbereich vorliegt!
leider trifft dieser Umkehrschluß zu !!!!!!!
wenn irgendein Vielfaches der Periodenlänge der Abgetasteten Funktion sich
der Zeit (samples/samplerate) also (Anzahl_Werte/Abtastfrequenz) nähert
dann werden die Spektralinien hoch und schmal, wenn dies nicht zutrifft
werden die Spektallinien entsprechend der Nichtübereinstimmung flach und
Breit, die Energie liegt also in der Fläche!
Außerdem steht zu befürchten das sich die Beiden Probleme überlagern
(belagern) weil sie nähmlich unter den selben Bedingungen auftreten.
Möglicherweise handelt es sich hierbei nicht um eine Überlagerung zweier
Probleme sondern lediglich um ein von Zeitbereich in den Frequenzbereich
transformiertes Problem, verwandt sind bei Probleme auf jeden Fall!
Es gibt zwei möglich Lösungen dieses Ärgerlichen Problems,
1. man benutzt eine geeignete Fensterfunktion windows(), daß heißt:
man multipliziert das Signal im Zeitbereich mit einer sanft Anklingenden
und sanft ausklingenden Funktion (sin,gauss...) so das zum Beispiel
eine Abgetastete Funktion nicht schlagartig anfängt und abbricht
sondern langsam ein und wieder ausgeblendet wird.
2. man nimmt das abzutastende Signal mehrfach mit unterschiedlichen
Abtastraten auf und überlagert die entstehenden Betragsspektren!
dabei empfehle ich instinktiv das Quadratische Mittel zur Überlagerung
heranzuziehen. Natürlich müssen die Frequenzen umgerechnet werden,
und genau da liegt ein weiterer gefährlicher Fallstrick, denn
die Mega-Blaster 16 (Mozart) Soundkarte emuliert ja nur eine
Soundblaster Pro Soundkarte, das heisst das sich aufgrund der
unterschiedlichen Hardware und des Begrenzten Wertebereich (8.Bit)
der Timerconstante kleine Unterschiede betreffend die tatsächliche
Abtastfreqenz ergeben. Die Folgen für ein Überlagerung wären Katast